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1 Abstract 

A detailed Monte-Carlo code has been developed from basic principles that simulates 
almost all of the basic photon and charged particle interactions. The code is used to 
derive the response functions of a high energy photon detector to incident beams of 
photons of various energies. The detector response matrice(DRM)s are calculated using 
this code. Deconvolution of an artificially generated spectrum is presented. 
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2 Introduction 

The objective of the present endeavour is to develop a simple and user-friendly detector 
simulation code. An output pulse height spectrum is generated by the present code that 
may be compared directly with a measured pulse height spectrum which is the observable 
distribution. 

The aim here is not to try to develop a simulation system as detailed and sophisticated 
as say, the EGSnrc or GEANT4 software systems but one that has all the basic physical 
processes incorporated within it and one which is, at the same time, sufficiently simple 
so that it may be used with ease to compare results of experiments where not very 
elaborate and complicated calibration systems are available- hence the calibration data 
also has certain inherent uncertainties and limitations due to effects (of say, surrounding 
materials etc.) that are neither well understood nor properly taken care of. In other 
words the present simulation code will be very useful in applications where extremely 
accurate results are not needed but a few percent accuracy will suffice. 

The approach taken here is to develop the entire detector simulation code starting 
from the very basic principles, generate artificial photon input energy spectra, let these 
interact with the detector to give rise to output pulse height spectra. Using the detector 
calibration data, the pulse height spectra may be converted to equivalent output energy 
spectra. Since the DRMs are calculated using the same program, the output energy spec- 
tra can be deconvolved using these DRMs and the deconvolved spectra can be compared 
with the original input photon spectra. This procedure establishes, to a large extent, the 
validity and accuracy of the entire simulation code. 

The simulation code is written in FORTRAN 77. The GNU compiler g77 is used to 
compile the code. It runs under Redhat Linux 9.0 operating system. A typical run to 
simulate a few thousand events for input photons say, 661 keV takes a few seconds in 
a PENTIUM IV system. This gives the energy deposition spectrum. In order to get 
the pulse height spectrum one has to use the photo-multiplier dynode multiplication 
simulation procedure. This takes rather long time, typically more than two hours for 
a few thosand events of 661 keV photons. The present code is nearly 3000 lines long 
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(excluding the data statements and the matrix inversion procedures). There are plans 
to rewrite the entire code in FORTRAN 90/95. 

3 The Detector Geometry and Co-ordinate Systems 

In the present simulations the detector geometry assumed is a cylindrical one. The center 
point of the detector is chosen to be the origin (O) of the global coordinate (right-handed 
rectangular Cartesian) system. The view axis of the detector is chosen to be the positive 
Z-axis of this co-ordinate system. 




z 




Figure 1: The global (OXYZ) and local (P'X'Y'Z') co-ordinate systems in the cylindrical 
geometry. 



3.1 Co-ordinate Transformations: 

Let a photon/particle start from a point P (x,y,z) and let its velocity (momentum) vector 
has a polar angle of 9 and azimuth of 4> with respect to the global co-ordinate system. If an 
interaction takes place at a geometrical distance dint{= dgml p) where dgm is the distance 
in gmcm~'^, p being the density of the medium, the co-ordinates of the interaction point 
P'(x',y',z') with respect to the global frame may be written as (Fig.l) 

X = x + dintsin6cos(f) (1) 

y =y + dintsinOsincf) (2) 

and 

Z = Z + dintCOsO (3) 

In the event of an interaction, let a product particle/photon be emitted at a polar angle 
9 and azimuth <j) where the Z -axis of the local co-ordinate system is along the velocity 
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(momentum) vector of the product particle/photon. The values of 9 and (f> in the local 
co-ordinate system are transformed back to the global system (O, $) as follows: 

sinQcos^ \ / cosOcosff) —sincj) sinOcoscj) \ I sin9 coscj) \ 
sinOsin^ = I cos9 sincj) coscj) sm9 sincj) I sin9' sincjl (4) 
cosQ ) \ -sin9 cos9 ) \ cos9' J 

Procedures have also been developed for rectangular geometry. Trial runs with these 
procedures gave encouraging results. 



4 Cross Sections 

Cross sections for all the different types of interactions of photons and charged particles 
are necessary to simulate these processes within the detector material. Since the compo- 
sition of the detector material is unknown a priori, cross section tables for all elements 
should be available to the program. At present the atomic cross section tables for only 
the following elements are available, viz. Hydrogen, Carbon, Sodium, Argon, Iodine, 
Xenon and Cesium. Near the absorption edges the cross sections are calculated for many 
closely spaced energy values. The XCOM program developed by Berger and Hubbell [1] 
is used for this purpose. This program is used to calculate the photo-electric, coherent, 
incoherent and pair production cross-sections. 

In the XCOM program, the coherent (Rayleigh) scattering cross sections are calculated 
using a combination of the Thompson formula and relativistic Hartree- Fock atomic form 
factors. 

In XCOM, the photo-electric cross sections are obtained following the method of 
Scofield [2] upto 1.5 MeV of input photon energy. At higher energies a semi-empirical 
formula [3] and the asymptotic high energy limit calculated by Pratt [4] is used. 

In XCOM, the incoherent (Compton) scattering cross sections are obtained using a 
combination of the Klein-Nishina formula and non-relativistic Hartree- Fock incoherent 
scattering functions. 

For pair production, XCOM uses combinations of formulas from Bethe-Heitler the- 
ory [5] alongwith other theoretical models to take into account screening. Coulomb and 
radiative corrections. Please refer to ref.l for the details of the XCOM program. Bremm- 
strahlung cross section is calculated separately. 

The total cross section for electron bremmstrahlung is taken as 

Trad = 4ao[to(183Z-l/3 _^ ^/^g] 

where cto = oiZ^r^, a is the fine structure constant, Z is the atomic no. of the material. 
Tg is the classical electron radius. Photon and charged particle interaction cross section 
tables are prepared for a wide range of energies, viz. 1 keV to 1000 MeV. 

Given the cross sections for the elements, the cross section for a composite material 
(such as Csl) for any process may be calculated as 

cr = T,niai (6) 

where n, is the number of atoms per unit volume of the ith element constituting the 
material. 

4.1 The Radiation Length 

The radiation length Xq is calculated as follows: 

1 /Xo = AaN/AZ{Z + l)r^/n(183Z- (7) 
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Here N is the Avogadro's number and A is the mass number of the clement. 
For a composite material the radiation length is given by the formula 

l/Xo=Pl/X,+P2/X2 (8) 

where pi, P2---, are the fractional weights of the various components and Xi, X2-.., are 
the corresponding radiation lengths. 

5 Simulations of Different Processes 

5.1 Photon Interactions 

For photons four different types of interactions are considered. These are: (i) Coherent 
scattering, (ii) Photo-electric absorption, (iii) Incoherent scattering and (iv) Pair pro- 
duction. Atomic relaxations are considered in the cases of Photo-electric absorption and 
Compton scattering. 

The geometrical distance of an interaction (d) from the point of incidence/previous 
interaction of the high energy photon is calculated from the equation 

d=-XhiR (9) 

where R is a uniform random number between and 1 . All the uniform random numbers 
between and 1 are generated using the fortran RAND() function. A is the mean free 
path (reciprocal of the total photon interaction cross section of the detection medium at 
the given energy). 

l/X = acoh. + <^p.e. + (^inc. + (^pair (10) 

where Uoo/i., <Tp.e., <^inc. and cTpair are respectively the coherent, photo-electric, incoherent 
and pair production cross sections. The value of the pair production cross section is, of 
course, equal to zero below a photon energy of 1.022 MeV. The type {K) of the interaction 
(one of the four types mentioned above) is determined using the relative (fractional) 
probability (cross section) for the specific process using a Russian Roulette, i.e. 

Pk = (^k/^(^k (11) 

cTfe and Pk being respectively the cross section and the relative probability for the process 
K. 

In the following we present the details of the algorithms used to simulate each of the 
above mentioned interactions. 

5.2 Photo-electric Absorption: 

In the photo-electric interaction with the atom, the incoming photon disappears and an 
electron is emitted whose energy equals the photon energy minus the binding energy of 
the electron in the relevant atomic shell. 

E = Ep-<j) (12) 

where the symbols have their usual meanings. The atom then undergoes relaxation 
through either fluorescent photon emission or non-radiative transition (described later). 
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5.3 Incoherent (Compton) Scattering: 

To simulate the incoherent scattering process first the atomic species is selected using 
a Russian Roulette (based on the relative incoherent scattering cross sections of the 
different constituent atoms of the detector material). Next, the relevant atomic shell is 
selected using the occupation numbers of the different atomic shells, ZijZ. The scattering 
angle (polar) Q of the photon in the local co-ordinates (the Z axis of the local co-ordinate 
system points along the momentum vector of the interacting photon) is generated using 
the Klein- Nishina formula for the differential scattering cross section 

da/rffi = Zrl{\l\ + a(l - x)f{\ + + c?{\ - xf l{\ + x''){\ + a(l - a;)]) (13) 

Here x = cos6, 6 being the polar angle of the scattered photon, a = 1/137 is the fine 

structure constant and is the classical electron radius. The von Neumann rejection 
technique is used for this purpose. Binding energy effects are taken into account by 
considering only those energy transfers to the electron that are larger than its binding 
energy. Dopplcr broadening effects are neglected. 

The energy of the scattered photon hi/ is calculated using the formula 

hu = hiy/[l + {hv/moC^){l - cosO)] (14) 

v being the frequency of the incident photon, h, mo and c are respectively the Planck's 
constant, the electron rest mass and the velocity of light. The azimuth <j) is uniformly 
generated between and 2-k. 

The polar scattering angle 6^ (in the local system) of the recoil electron is calculated 

as 

6»e = tan"^ (1/(1 -I- a)tan{QMp)) (15) 

where Bp is the scattering angle of the photon and a = hu/mo(?- The azimuth cj) of the 
recoil electron is equal to {w — (j)p) while its energy is given by 

Ee = E-Ep (16) 

5.4 Coherent Scattering: 

The angular distribution ducoh/d^ for coherent (Rayleigh) scattering is given approx- 
imately, by the distribution function dar/dO. for classical Thomson scattering by an 
electron, 

daT/dVl = rl/2{l + cos^e) (17) 

5.5 Electron-Positron Pair-production 

The screening parameter 7 is defined as 

7 = im{moC^/E)[l/v{l - v)]Z-^'^ (18) 

where 

v = {E' +moC^)/E (19) 

Here E is the energy of the photon that undergoes pair production into a positron having 
total energy E and an electron having total energy E — E . nio is the electron rest 
mass, c the velocity of light and Z is the atomic number of the target nucleus. The 
differential cross section for pair production ^pair{E, E ) depends on the value of the 
screening parameter. The theoretical expression for ^pair(E,E') given by [6] as 

%^ir{E,E')dE' =4aN{ZyA)rl{dE' /E)G{E,v) (20) 
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For different ranges of the value of the screening parameter, the following expressions are 
used. No screening (7 >> 1): 

G{E, v) = [v"^ + (1 - vf + {2/3)v{l - v)][ln{2E/moC^)v{l - v) - 1/2] (21) 

Complete screening (7 = 0): 

G{E, V) = [v^ + (1 - vf + (2/3)t;(l - v)]ln{183Z-^/^) - (l/9)u(l - v) (22) 

Intermediate cases (0 < 7 < 2): 

G{E,v) = [v' + (1 - vnh{j)/4-{l/3)lnZ] (23) 

(2 < 7 < 15): 

G{E,v) = + (1 - vf + (2/3)t;(l - v)][ln{2E/moc'^)v{l - v) - 1/2 - 0(7)] (24) 

The functions /i(7), /2(7), and 0(7) are the same that enter in the expressions for the 
radiation probabilities mentioned in subsection 7.1. For very small energy, k very near 
to 2moC^, Hough's approximate formula 

^E+)dE+ = 8/3$o(fc - 2moC^lkfzdE+ (25) 

is used, where 

z = 2y/x{l - x) (26) 

and x is given by 

x = E+ -■moC^/k-2moC^ (27) 

The angle of the positron (electron) with respect to the momentum vector of the high 
energy photon (Z' axis of the local system) is given by 

9 = moC^/E± (28) 
5.6 Muon Pair Production: 

The formulae given in the previous subsection may be used to simulate /i+ — /i~ pairs 
although the threshold for this process is equal to 206 times the threshold for the corre- 
sponding process for e+ — e~ pair, i.e. 210.532 MeV. 

6 Atomic Relaxations: 

In the case of photo-electric absorption and also in the case of incoherent scattering, a 

vacancy is created in the electronic shell with which the photon interacts. This vacancy 
is filled by electrons from higher energy shells and the excitation energy is released in the 
form of either (i) one or more fluorescent photons, (ii) an Auger electron, or (iii) a Coster- 
Kronig electron. The relative probabilities for these three types of relaxation processes 
are different for different electronic shells. In the present work the Coster-Kronig process 
is neglected. 
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6.1 Fluorescence 



The energy of the fluorcsnce photon is equal to the difference in the binding energies 
of the initial and final states of the atom. The fluorescent X-ray photon is emitted 
isotropically. Its polar angle {6) and azimuth (</>) in the local system are calculated as 
follows: 

e = cos-\2Ri - 1) (29) 

and 

(j) = 27ri?2 (30) 

Ri and R2 being two different random numbers that are distributed uniformly between 
and 1. The values of 9 and (p are transformed back to the global co-ordinate system 
using eqns.(l) and (2). 

The fluorescent photon may, either (i) escape the detection medium (the location of 

its interaction being outside the boundary of the detector and in that case its energy is 
treated as lost resulting in an energy deposit that is equal to the input photon energy 
minus the fluorescent photon energy. This kind of events give rise to the characteristic 
X-ray escape peak), or, (ii) it may interact photo-electrically with a higher shell in which 
case either its total energy or a part of its energy (in the case of escape of a higher shell 
fluorescent photon) is deposited. 

6.2 Auger Electron Emission 

The Auger electron emission (non-radiative electron emission) is considered only for the 
K and L shells. The energy(ies) of the Auger electron(s) is completely absorbed within 
the detection medium and contribute to the energy deposited. 

7 Charged Particle Interactions: 

For electrons, bremmstrahlung and for positrons, brommstrahlung and pair annihilation 
processes are considered. For electrons and positrons multiple Coulomb scattering is 
considered in an approximate manner. 

7.1 Bremmstrahlung 

The Bremmstrahlung process is essentially similar to the pair production phenomenon. 

The differential cross section (probability) for the radiation process, ^rad{E,E ) de- 
pends on the value of the screening parameter which is deflned as 

7 = 100(TOocVC/)b/(l - v)]Z-'^/^ (31) 

Here, 

U = E + ruoc^ (32) 

and 

V = e'/U (33) 
The expression for $rad(-S> E ) is given by ref.6 as 

^rad{E,E')dE' =4a{N/A)Z^rlidE'/E')F{U,v) (34) 

For different ranges of the value of 7, the differential cross section is given by the 
following formulas: No screening (7 >> 1 ): 

F{U, v) = [l + {l- vf - (2/3)(l - v)][ln{2U/moc''){l - v)/v - 1/2] (35) 
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Complete screening (7 = ): 

F{U, v) = [l + {l- vf - (2/3)(l - v)]ln{imZ-'^'^) + 1/9(1 - v) (36) 
Intermediate cases, (0 < 7 < 2 ): 

F(C/, = [1 + (1 - vf\[h{^)l^ - (l/3)ZnZ] - (2/3) (1 - v)[h{l)l^ - (l/3)^nZ] (37) 
(2 < 7 < 15 ): 

F{U,v) = [1 + (1 - vf - (2/3)(l - v)][ln{2U/moC^){l - v)/v - 1/2 - c{^fl^] (38) 

The average angle of the emitted photon with respect to the momentum vector of the 
high energy electron (Z axis of the local system) is given by 

^0 = moC^/Eo (39) 

7.2 e^e~ Annihilation 

Only electron-positron pair annihilation at rest is presently considered. The energy of 
each annihilation photon is equal to 0.511 MeV. The direction of the second photon is 
taken to be exactly opposite to the direction of the first photon which is sampled from a 
distribution that is isotropically distributed. 

7.3 Knock-On {5 ray) production 

High energy electrons often produce electrons of energies that are comparable to the 
energies of the primary electrons. This process is simulated using the formula [7] 

P{E)dE = WdEjE^ (40) 

where P{E)dE is the probability of an electron receiving the energy E. is a constant 
that depends on Z , A\ the atomic number and the mass number of the medium in 
addition to the velocity /? of the primary electron. 

7.4 Multiple Coulomb Scattering 

The well-known Moliere theory of Coulomb scattering is used in order to take into account 
the scattering of light charged particles (electrons and positrons) in the Coulomb fields 
of both atomic nuclei and atomic electrons. 

For this the energy domain is divided into two intervals. The first interval includes 
charged particles whose energies lie between 1 kcV and 1 MeV. The second interval 
includes the energy region between 1 MeV and 1 GeV. These two energy intervals are 
again subdivided each into 30 logarithmic energy bins (following the treatment used by 
Vatcha [8]). For charged particles having energies above 1 GeV, only a gaussian having 
a width equal to [Esj Efi)s/(t), (where E is the energy of the particle, (3 = vIc,y being 
its velocity and t = x/Xq\s the thickness of the slab of material in terms of its radiation 
length) is considered. The constant Eg = 21 MeV. 

Bethe [9] has expressed Moliere's theory of multiple Coulomb scattering in a form 
which is easy to simulate. Using the first two terms of his equation () together with a 
correction for solid angle the frequency function for scattered angle is 

f{@)d@ = (sme/e)V2[/(o)(0) + l/B/(i)(0)]0d</, (41) 
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Here <j) is the reduced angle of which is defined as 



= e/XcSV2 (42) 

where 

xl = Q.l^7Z{Z+l)t/A{iwf (43) 

where pv is in MeV, t is in grams per square centimeter, and A is the atomic weight of 
the material. The parameter B is defined by the transcendental equation 

B-lnB = B (44) 

where 

B = 2ln{xc/Xa) + 1 - 2C (45) 

where C = 0.577... is Euler's constant and Xa is an angle which depends on the screening 
of the atomic field. It is defined as 

xl = {X/5MlaoZ-^/Y{l.U + 3.76a^) (46) 

with 

a = 2TrzZ^/hv (47) 
The functions f'-°\(j)) and f''^\(j)) are defined as 

/(o)(<^) =2e-^' (48) 

and 

= f^'Hm' - i)m^'') - Inicf") - 2(1 - /(°)(<^)) (49) 
where Ei{x) is defined as follows: 

/inf f^x 
e-*/tdt= / e*/tdt (50) 
-X t/ — inf 

For a given electron energy a slab thickness (proportional to the energy of the electron) 
of the material is chosen. The percentage probability of scattering is selected using a 
uniform random number between and 1. The angle of scattering is selected from the 
table using reverse interpolation [10]. 

7.5 lonisation Loss and Particle Ranges 

For the ionisation loss the expression given by Bethe and Bloch in ref.7 is used. 

dE/dx = {l/p)dE/dX (51) 
Here x = Xp is in gmscm~^ while X is in cms. 

dE/dX = -K{Z/A){p/l3'^)[ln2moC^I3'^EM/l'^{l - (3'^) - 20^] (52) 

where 

K = 2TrNz'^e^/moC^ (53) 

Here N is the Avogadro number, nio and e are mass of the electron and its charge, Z, 
A and p are the atomic number and the mass number and the density of the medium, 
respectively, and / is its effective ionization potential; z is the charge and /3 the velocity 
(in units of c, the speed of light) of the incident particle. 



9 



Em represents the maximum possible energy transfer in an interaction and is given 

by 

Em = 2moc2/3V(l - (54) 

For energies of electrons {6 rays) an approximate formula for the practical range, in 
gcm~'^ is used (ref.7) 

Rp = 0.71£;i-^2 (55) 
where the energy of the electron (E) is expressed in MeV. 

7.6 Emission of Scintillation Photons 

The number of scintillation photons (integrated over the spectral band) produced in the 
scintillator is obtained simply by dividing the energy dcposicd by a constant that gives 
the average energy rcqired to produce a scintillation photon. The value of this constant 
is 0.02564 eV for Csl [11]. Since the number of scintillation photons is always an integer, 
0.5 is added to the calculated value and the result is truncated. 

7.7 Emission of Cerenkov Photons 

In case the charged particles produced in any of the interactions happen to traverse 
regions that consist of transparent media (such as air, water, plastics etc.) and the 
velocity of the particles exceed the velocity of light in those media, Cerenkov photons 
are emitted. 

The angle of emission 9c of the Cerenkov photons is given by the equation 

COS0C = l/n/3 (56) 

where [3 is equal to v/c, v being the velocity of the particle, while n represents the 
refractive index of the local medium. The number of Cerenkov photons dN emitted in 
the track length element dl in the wavelength interval AA is given by the equation 

dN = 2Tradl{l-l/n'^l3^)AX (57) 

The spectral distribution of Cerenkov photons is sampled using the probability distribu- 
tion 

P{X)dX = {1.0/X'^)dX (58) 



7.8 Chcirge Multiplication at the Photomultiplier: 

The scintillation photons incident on the photo-cathode (Ni) of the photo multiplier 
(PMT) are converted into a photo-current (number of photo-electrons, A^e) 

Ne=<v> Ni (59) 

where < r] >, the average spectral efficiency of the photo-cathode is taken to be equal 
to 0.2. < 77 > depends on the spectral response of the particular photocathode type. 
Strictly speaking, 

iVe = 1 N,iX)vW/ 1 (60) 

where Ni{X) is the number of scintillation photons having wavelength A and t]{X) is the 
corresponding spectral efficiency. 

The photo-electrons are multiplied by the successive dynodes (8 to 10 in number) 
through secondary emission process. 
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The energy resolution of the detector is determined largely by the statistical fluctu- 
ations in the minimum number of particles (electrons) in the dynode chain, i.e. the 
number of photo-electrons. The charge multiplication at the dynodes is simulated in the 

following manner. 

Using the knowledge of the total PMT high voltage and its distribution between 
different dynodes (usually all the dynode voltages are equal except the first one), from the 
PMT data book the average secondary electron multiplication factors are obtained. Since 
the electron multiplication is a Poisson process, in each succesive stage this multiplication 
is simulated using a previously calculated table of Poisson probabilities for obtaining an 
integral number of secondary electrons for a given average value. 

7.9 Folding the Radial Response in the case of scintillation de- 
tectors 

It has been observed that in scintillation detectors, the photomultiplier tube light collec- 
tion efficiency reduces to 85% near the; edge of the crystal for a given energy deposition, 
as compared to the same energy deposited at the center of the crystal. This radial radial 
response function (a quadratic fit to light collection vs. radial distance of the point of 
energy deposition) is used to obtain the amount of light incident at the photo-cathode 
of the PMT. 

7.10 Chcirge Integration and Pulse Height Spectrum 

The charge output at the anode of the PMT is integrated into a charge sensitive pre- 
amplifier (CSPA) having a charge gain equal to 0.025X10^^ Volts/Coulomb, thus pro- 
ducing a voltage pulse as output. This is amplified using a linear amplifier of gain equal 
to 5. 

7.11 Calibrating the Detector 

The detector is calibrated as follows: The response functions (the pulse height spectra) 
are simulated for two monoenergetic incident photons (say, 100 keV and 600 keV) using 
a large (say, 10 000) number of incident photons. The positions of the photo-peaks (full 
energy peaks) are determined accurately using very fine binnings. A straight line fit is 
obtained with the pulse height as a function of the input photon energy 

h = mE + b (61) 

h being the pulse height in volts and E being the photon energy in keV. The calciilatcd 
values of m and b are used to convert output pulse heights into corresponding energy 
values to form the output energy spectra. 

8 Generation of Artificial Photon Input Spectrum and 
Spectral Deconvolution: 

In order to check the correctness of the calculated DRMs it is necessary to generate 
artificial photon energy spectra, let them interact with the detector and obtain output 
pulse height spectra. These pulse height spectra should then be deconvolved back into 
photon energy spectra using the calculated DRMs. These deconvolved photon spectra 
should be compared with the original photon spectra using say, tests. This procedure 
would then establish the correctness of the entire detector simulation code. 
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8.0.1 Power- Law Photon Energy Spectrum 

A power-law (differential) photon energy spectrum having the form AE~'^ where A is the 
normalisation constant, E the energy of the photon and 7 is the spectral (differential) 
index, is generated using the inverse transform method as follows: 

E={KR + kiy^^^-^^ (62) 

where K = £^2^ ''^^ — ki and ki = ''K Ei and E2 are respectively the minimum 
and maximum limits of the spectral energy band, i? is a random number uniformly 
distrubuted between and 1. 

An input spectrum of photons (power law) is simulated. The corresponding output 
pulse height spectrum is converted to an equivalent energy spectrum (using the artificial 
calibration data). Let P be the vector that represents the output pulse height spectrum 
(rebinned) and R the detector response matrix. If the true input photon energy spectrum 
is denoted by the vector S, then 

Pi = SRijSj (63) 

where Pi is the detected counts in the ith energy channel, Sj the number of input photons 
in the jth energy channel and Rij is the ijth clement of the DRM. The true incident 
photon spectrum S is obtained by the matrix inversion procedure. 

S = R-^P (64) 

where R"-*- is the inverse of the detector response matrix (historically, the first such 
detector response matrix was calculated by J. H. Hubbell [12]). 



8.1 Simulation of Gamma Ray Burst (GRB) Spectrum 

The energy spectrum of Gamma Ray Burst sources are best described by the Band 
[13], [14] spectrum. A typical Band spectrum is given as 

N{E) = A{E/ 100)" exp{-E/Eo) (65) 

for 

E<={a- p)Eo (66) 

and 

N{E) = A{a - P){Eo/100)°'-^{E/100f (67) 

for 

E>{a- P)Eo 



9 Gamma Ray Applications 

In the case of gamma ray studies, the detector (e.g. a Csl scintillator) is usually shielded 

using a plastic scitillator that is used as a veto counter to remove the charged particle 
background. Usually there is a thin window (say, of Alumnium) in front of the scintil- 
lator. In addition a metal casing is used to cover the detector. Gamma rays interact in 
all of these materials and produce either coherently scattered gamma ray photons hav- 
ing energies equal to that of the primary, or, they may produce incoherently scattered 
photons of lower energies or they may also produce annihilation gamma rays. These 
secondary photons have finite chances of entering the actual detection volume (in this 
case the Csl crystal). Sometimes they may be scattered from another material and enter 
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the scitillator though the probabihties of such tertiary and higher order scatterings are 
small. 

These effects are taken care of in the following manner. Nested volumes of materials 
arc considered where the actual detection volume (in this case Csl) is designated as the 
volume of zeroth order. The next set of volumes (composition will in general be very 
different) that enclose the zeroth volume is called the volumes of first order and so on. 
When a primary gamma ray is incident on any of the volumes, its history is traced until 
either it misses or hits the zeroth volume. If it misses the zeroth volume it is rejected 
while in the other case it is considered as a valid incident photon of degraded energy. 
Development of this procedure is yet to be completed. 



10 Results 

In the following the results of simulations for a Cesium Iodide scintillation detector (3 
inch diameter and 0.5 inch thickness) are presented. 

10.1 Detector Efficiency 

The efficiency of a Cesium Iodide scintillation detector has been estimated from the 
simulations. These estimated values of the detector efficiency are plotted in Fig. 2. 

EFFICIENCY OF CESIUM IODIDE DETECTOR 



100 200 300 400 500 

Photon Eneigy (keV) 



Figure 2: The efficiency of the Csl detector (3 inch diameter and 0.5 inch thickness) is 
plotted as a function of energy. 



10.2 Energy Resolution of the Detector 

From the simulations the energy resolutions of a Cesium Iodide scintillation detector has 
been derived at various energy values. The energy resolution is obtained by dividing 
the full width at half maximum (FWHM) of the photo-peak (full energy peak) by the 
energy of the photo peak. In Fig. 3 the estimated energy resolutions are plotted against 
the energy of the incident high energy photon. 

10.3 Sample Response Functions 

In Fig. 4 a sample response function is shown. The energy of the incident photons are each 
equal to 70 keV. In Fig. 4 the peak near 20 keV is due to the escape of the characteristic 
X-ray ffuorescence photons. For Cesium and Iodine, the energies of these photons are 
respectively 35.98 and 33.17 keV. In Fig. 5 the response function of the same detector for 
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ENERGY RESOLUTION OF CESIUM IODIDE DETECTOR 



10 100 1000 

Log (Energy (keV)) 

Figure 3: The energy resolution of the Csl detector is plotted as a function of energy. 
The E~^l'^ dependence of the energy resolution is clearly evident in the figure. 



RESPONSE FUNCTION OF Csl DETECTOR FOR 70 keV INPUT PHOTONS 
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Figure 4: Simulated response function of the Csl scintillator to an input beam of 70 keV 
photons. The K X-ray escape peak is clearly seen. 
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525 keV input photons are shown. 



RESPONSE FUNCTION OF Csl DETECTOR FOR 525 keV INPUT PHOTONS 
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Figure 5: Simulated response function of the Csl scintihator to an input beam of 525 
keV photons. The Compton edge and the Compton continuum are clearly visible. 



10.4 Response Matrix and Spectral Deconvolution 

The response functions of the Csl detector are calculated for the following 16 different 
energy values, viz. 28, 44, 60, 76, 92, 120, 160, 200, 240, 280, 325, 375, 425, 475, 525 
and 575 keVs. These response functions are binned into the following energy intervals, 
viz. 12-30, 30-49, 49-68, 68-87, 87-106, 106-152, 152-199, 199-246, 246-293. 293-340, 340- 
398, 398-457, 457-516, 516-574 and 574-633 keV respectively. Thus, the response matrix 
(DRM) of the Csl detector is obtained. 

A power-law (differential index equal to -1.38) photon spectrum is generated using 
the algorithm described in section 7.0.1. A total of 9 936 photons having their energies 
between 20 keV and 600 keV are generated. In addition, 70 photons having energy equal 
to 511 keV each are generated. These photons are made incident on the Csl scintillator 
at an angle of 1 degree relative to the view axis. These photons interact within the 
scintillation crystal and give the output spectrum depicted in fig. 6. This output spectrum 
is obtained by binning the output energy spectrum using the energy intervals described 
in the previous paragraph. 
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Figure 6: Simulated output spectrum from the detector for an input power law photon 
spectrum (differential index equal to —1.38. At higher energies the counts decrease due 
to the gradually decreasing efficiency of photon detection. 
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The inverse of the DRM is calculated using well established procedures [15]. The 
result of deconvolution of the output energy spectrum of fig. 6 with this inverse of the 
DRM is shown in the following figure. 



SPECTRAL DECONVOLUTION 



Figure 7: The deconvolved energy spectrum (crosses) is compared with the input photon 
spectrum (plus signs). The two spectra agree quite well at low and medium energies but 
seem not to match very well at the high energy end. 



11 Results of Simulations of Other Interaction Pro- 
cesses 

Pair production at high energies is simulated using the formulae given in section 4.5. The 
results are shown in the figure below. In Fig. 8 the fractional energy (v) of the positron is 
plotted along the abscissa and the quantity E(j)pair {E, E ) is plotted along the ordinate. 
Here (l)pair{E,E ) = Xo^pair{E,E ) is the differential pair production probability per 
radiation length, ^pair {E, E ) is the differential pair production probability of production 
of a positron having energy E by an incident photon having energy E in IgmcmT^ 
thickness of the material. 



Differential Pair Production Probabilities 
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Figure 8: Simulated differential pair production probabilities at high energies. The values 
(inset) give the energies of the input photons in terms of the electron mass. 



The pair production at low energies are simulated using the Hough's formulae given 
in section 4.5. The results of one such simulation is shown in the following figure. 
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Differential Pair Production Probabilities 
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Figure 9: Simulated differential pair production probabilities at low energies. The value 
(inset) give the energy of the input photons in terms of the electron rest mass. 



The electron/positron brcmmstrahlung processes at high energies arc simulated using 
the formulae given in section 6.1. The results of these simulations are given in fig. 10. Here 
again the fractional energy (v) of the emitted photon is plotted along the abscissa and the 
quantity E (f>rad{E, E'), where (l>rad(E, E ) = Xo^j-adiE, E ), is the differential radiation 
probility per radiation length. $pair {E, E ) is the differential radiation probability of 
production of a photon having energy E by an incident electron having energy E in 
lgmcm~'^ thickness of the material. 
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Figure 10: Simulated differential radiation probabilities at high energies. The values 
(inset) give the energies of the input electrons in terms of the electron rest mass. 



Knock-on (delta ray) productions by a 10 MeV electron are simulated using the pro- 
cedure described in section 6.3. The results are shown in fig. 11. 

Ccrenkov photon emission by high energy electrons/positrons are simulated using the 
formulae given in section. 6. 7. The simulated wavelength spectrum of Cerenkov photons 
in the wavelength band of 300nm to 650nm is shown in the following figure. 

12 Simulations of Photon Energy Spectra 

Different sources emit different types of high energy photon spectra. It is necessary, 
therefore, to be able to simulate different types of photon energy spectra. In the following 
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Delta Flay Production from a 10 MeV electron 




Figure 11: Simulated probability distribution for the production of delta rays (knock-on 
electrons) produced by an incident electron of energy 10 MeV. 



CERENKOV PHOTON EMISSION SPECTRUM 
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Figure 12: Simulated wavelength spectrum of cerenkov photons emitted by an elec- 
tron/positron within the wavelength range of 300nm to 650nm. 
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two different types of spectra, viz. (a) a power-law energy spectrum (large number of 
photons- the total being equal to 9 million) and (b) a smoothly-broken power-law (Band 
Spectrum) spectrum are shown. 




Figure 13: Simulated power-law energy spectrum of photons having a differential energy 
index of —1.38. 



Gamma Ray Burst Speclrui 




Figure 14: Simulated energy spectrum (Band spectrum) of a hypothetical gamma ray 
burst. The parameter values are: (i) A = 1.0S4, (ii) a = —1.0, (iii) /? = —2.0 and (iv) 
Eq = 150.0keV. 



13 Discussion 

It is possible to incorporate any detector geometry (say, rectangular geometry) in the 
present simulation code and calculate the response of that detector. Also, the response 

functions of other types of detectors, for example, gas (Argon or Xenon) filled multi-wire 
proportional chambers (MWPCs) such as those used in X-ray astronomy experiments 
may be obtained. 

In the present simulations the differential cross sections used for the coherent scatter- 
ing and the incoherent scattering processes are respectively the classical Thomson cross 
section and the Klein-Nishina formula. Strictly speaking, in the case of coherent scatter- 
ing the square of the atomic form factor F{q, Z) should be multiplied with the Thomson 
cross section while in the case of incoherent scattering the incoherent scattering function 
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S{x, Z) should be multiplied with the Klein-Nishina differential cross section. These 
changes will be introduced in a later version of the code. 

All the different modes of atomic relaxations are yet to be incorporated in the present 
code. Pair annihilation in flight and triplet production are also to be introduced. 

There are some small irregularities and discontinuities in the plots of the differential 
probabilities for pair production and bremmstrahlung (fig. 5 and fig. 7). Presumably these 
have resulted due to the fact that somewhat coarse interpolations of the three screening 
functions (/i(7), f2{l), and 0(7)) have been used in the present simulations. 

There is minor disagreement of the deconvolved spectrum with the input photon spec- 
trum (Fig. 5) at the high energy end. This is due to the fact that the number of events in 
this energy range are small and the response functions at the higher energies are poorly 
determined (due to the insufRency of the number of events generated to calculate the 
response functions). 

14 Conclusions 

A high energy photon detector simulation code developed from the first principles has 
been described. The results obtained using this code are quite encouraging. This code 
gives an output pulse height spectrum, instead of an energy loss spectrum. This simula- 
tion code is being developed further to incorporate all the possible interaction types of 
photons and particles. It should be possible to extend this code to simulate extensive 
air showers (EAS) in the atmosphere produced by very high energy (VHE) and ultra 
high energy (UHE) gamma rays and electrons and also to study the cerenkov radiation 
produced in these cascades. 
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